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The Smoluchowsky equation for a system of interacting Brownian particles in a temperature 
gradient is derived from the Kramers equation by means of a multiple time-scale method. The 
interparticle interactions are assumed to be represented by a mean-field description. We present 
numerical results that compare well with the theoretical prediction together with an extensive dis- 
cussion on the prescription of the Langevin equation in overdamped systems. 
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I. INTRODUCTION 

The erratic motion of tiny particles suspended in a fluid of lighter particles is a known phenomenon, going under 
the name of Brownian motion, of fundamental interest in Statistical Mechanicsi*^. A classical method to study 
this physical process is the Langevin equation, where the influence of the solvent particles on the heavier particles is 
modeled by two force terms. One term represents a viscous force, with friction coefficient 7, the other a stochastic force, 
whose intensity depends on the temperature of the solvent. Alternatively, the description of the dynamics of a system 
of Brownian particles in an external field may rest upon the Fokker-Planck equation (FPE) for the probability density 
of the particles. The two levels of description, namely, a) the Klein-Kramers equation^, which considers the evolution 
of the probability density of position and velocity, and b) the Smoluchowski equation^, for the probability distribution 
of position only, have deserved a great deal of attention. The passage from the first description to the second one 
is well understood for the case of ideal non-interacting particles immersed in a heat bath at constant temperature. 
In fact, as shown by Wilemski and Titulaer 6 , the derivation of the Smoluchowski equation from the Kramers-Klein 
equation, can be performed by means of a perturbation expansion in terms of the inverse of the friction coefficient, 
7. On the other side, one can derive the Smoluchowski equation more directly after taking the overdamped limit 
7 — > 00 in the Langevin equation and neglecting the inertial acceleration term. Physically speaking, this procedure 
corresponds to the fact that the velocities of the particles thermalize rapidly, so that only the probability distribution 
of positions of the particles matters. 

Thermal gradients, besides causing a heat flow, can induce a mass flow in systems such as colloids. The phe- 
nomenon is termed thermodiffusion, thermophoresis or Ludwig-Soret e ffectp£i2*2*i2iiiii2iiii and is relevant because it 
represents a tool to manipulate and concentrate molecules in solution. Thermophoresis, in fact, typically enhances 
the concentration in colder regions. Recently, it has been employed to form two dimensional colloidal crystals using 
1\xm polystyrene beads. Such a method leads to applications in microfiltration, particle accumulation and molecular 
detection on surfaces^. Nevertheless, in spite of its technological interest, not many studies provide a derivation of 
the governing equation for the concentration field in the presence of a temperature gradient externally imposed. A 
problem arises when one considers the 7 — > 00 limit: the resulting stochastic equation contains a multiplicative noise 
term (the temperature is space-dependent). This in turn means that the associated FPE depends on the particular 
prescription used to derive it, i.e. one encounters the so called Ito-Stratonovich dilemma. To solve it, Matsuo and 
Sasaiii recently obtained the Smoluchowski equation from a perturbation expansion in powers of the inverse of 7. 
Two important messages arise from their work: first, the underdamped dynamics is free from the Ito-Stratonovich 
dilemma, and, second, their perturbation scheme is unbounded in time and one needs a renormalization treatment to 
make it convergent. 

The main objective of this work is to present a systematic and convergent derivation of the Smoluchowski equation 
for a system of interacting Brownian particles in a temperature gradient. The interaction is of mean-field type, so 
that our starting point, the Kramers equation, can be written down in close form for the one-particle distribution 
function. We employ a multiple time-scale method which has been designed to deal with non-uniformities in systems 
with more than a time scale. The approach is quite appropriate for our objectives, since in the large friction limit 
there is a very clear time-scale separation in the system: in a very short period the velocities of the particles relax 
to their thermal equilibrium value, whereas on a much longer time scale also the positions of the particles assume 
configurations compatible with their steady state distribution. In this regard, the multiple time-scale method has been 
used to perform the passage from the Kramers to the Smoluchowski equation for an ideal gas at constant temperature 
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as an expansion in powers of the small parameter 7 -1 , providing a uniformly valid resul1pi£iii. More recently, the case 
of interacting hard-core particles in isothermal systems has also been studied within this frameworki*. 

Another important part of this work is the discussion about the Ito versus Stratonovich description of the system. 
This is based on numerical results of the Langevin dynamics. 

The plan of this paper is the following: in section [H] we present the full Langevin equation, its overdamped limit, 
and the corresponding Kramers and Smoluchowski descriptions of the system. In section IIIII the derivation of the 
Smoluchowski equation with the multiple time-scale method is presented. In section IIVI we present the numerical 
results, with some discussion on the different Ito-Stratonovich prescriptions, and in section[V]we give our conclusions. 



II. MICROSCOPIC MODEL. SMOLUCHOWSKI AND KRAMERS EQUATIONS. 

Let us consider a system of N interacting Brownian particles of mass m in an external force field. For the sake of 
simplicity, let us assume only one spatial dimension. The system is also under the influence of an external non-uniform 
temperature field so that the dynamics of the particles is given by: 



m llt = -mjVi + fext(Xi) + ^2 fint(Xi - Xj) + y/T(Xi)^i(t), (1) 

for i = 1, ...,N, and where 7 is the friction coefficient, T is a non-uniform heat-bath temperature, f ex t is an external 
force, fi n t is the interaction force between pair of particles, and is a Gaussian white noise with zero average and 
correlation 

(&(*)&(*)) = 2"/mk B S l3 6{t - s). (2) 

(•) indicates the average over a statistical ensemble of realizations, and ks is the Boltzmann constant. To study the 
dynamics of the system we consider the single particle probability distribution function, P(xi,Vi,t), of the position 
and velocity variables. Its evolution is given by the Fokker-Planck equation^ also known as Kramers equation: 



d_ 

dt 



p(x,v,ty 



fext(x) + f m (x) 



d v 



P(x, v, t) = jd v vP(x, v, t) + kBT ^ d v P(x, v, t) 



Here f m (x) is the molecular field which takes the form 



(3) 



(4) 



being U(\x — x'\) the pair potential energy (i.e. fint( x i ~ x j) = 8 g* . ); and p(x) = f dvP(x,v,t). It is very 
important to realize under what conditions on the interparticle interactions one can write down eq. J3Jl (see**). Note 
that with no approximations the Kramers-FPE for interacting particles for the one-particle distribution function, 
P(x,v,t), depends on the two-particle distribution P^ix, v, x', u', t). The main hypothesis done in Eq. © is that 
whenever U[x — x') is a smooth function of particle separation one may assume a mean- field approximation for the 
two-particle distribution function P<z{x, v, x',v', t) » P(x, v, t)P(x' , v', t). In the following we shall define / = fext+fm- 
Since the temperature T(x) in Eq. @ multiplies a derivative with respect to v, one realizes that the Kramers 
equation is free from the so called Ito-Stratonovich dilemma. On the contrary, when one considers the overdamped 
limit, 7 — ► 00, in Eq. one obtains the following Langevin equation: 



dxi 
~~dt 



f( Xi ) , 



TO7 



(5) 



which leads to different equations for the distribution of the spatial position (the so-called Smoluchowski equation) 
depending on the prescription adopted to integrate the stochastic evolution. In fact, the corresponding FPE for Ito 
and Stratonovich conventions, the two rules usually applied, are respectively: 



P(x,t) = -^d x (fP) + ^d*(TP), 



P(x,t) 



III 



dx(fP) 



d x {T^d x {T^P)). 



(0) 
(7) 
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The aim of this work is to derive a Smoluchowski equation for the system using the multiple time-scale method, 
via an expansion in the inverse of the friction coefficient. Our starting point is the Kramers equation, Eq. @. At 
this stage it is convenient to introduce the following dimensionless variables: 

Vt v x I 

r = t-f, V=—, X = ~, r = 7 — , 8 

lf(x) 

F(X) = P(X, V, t) = lv T P(x, v, t), (9) 

mv T 

being vt = -^fcsTo/m the thermal velocity at the reference temperature To (which can be taken equal to 1), and I is 
a typical length scale of the system, for example the effective radius of the particles. The non-dimensional Kramers 
equation for P(X, V,r) denoting, f(X) = T(x), is: 

dP(x v,t) = Tdv , v pj + T fr X )dlP - vd x P - Fd v p. (io) 

OT 

In the following we shall skip the tildes in our notation. The (non-dimensional) Kramers equation, Eq. (|10|l . is our 
definitive starting point and in the next section we proceed to present the set of eigenfunctions that will be used in 
the subsequent section when we make use of the multiple time-scale approach. 

III. METHOD OF SOLUTION 

Standard solutions of Eq. go through an expansion in a basis set of eigenfunctions. First, it is convenient to 
separate the spatial dependence from the velocity dependence of the probability distribution. We rewrite Eq. (|10|l as: 



TL fp P(X,V,t) 

where we introduced the Fokker-Planck operator 



[ — + V—+F— 
8t dX dV 



P(X,V,t), (11) 
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dV 



d 



dV 



Lfp P{X,V,t) = — T(X)— + V P{X,V,t). (12) 



We will expand P in terms of the eigenfunctions H„(X, V) of Lpp^ which are related to the Hermite polynomials: 



L FP H„(X, V) = -vH v (X,V) 



(13) 



with 



Tixyi 2 d v V 2 



(14) 



We now express the solutions to Eq. (|11|) as a linear combination of products of the type <P^(X,t)H^(X,V) and 
after some algebra obtain: 



v=0 



(15) 



T'(X) 



b v {X,r) F{X) 



MX,t) H„ +1 {X,V) + ^T(X)v- 



,(X,t) 



zvnx) 



ox T{xy" y ' '\ " +1V ' ' v v ' ox 

MX, t)[H v+3 (X, V) + 2(1/ + 1)H V+1 (X, V) + vH v _ x {X, V)]}= 



(^o-l)^-i(^, V) 



where the prime denotes derivation with respect to the argument. Equating the coefficients of the same basis functions, 
H v , we obtain an infinite hierarchy of equations relating the various moments 4> V {X 1 r). Such a hierarchy can be 
truncated by the multiple scale method, as shown hereafter. 
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A. Multiple time-scale analysis 



In the multiple time-scale analysis one determines the temporal evolution of the distribution function P(X, V, r) in 
the regime << 1 by means of a perturbative method. 

In order to construct the solution one replaces the single physical time scale, r, by a series of auxiliary time scales 
(to, t\, .., r n ) which are related to the original variable by the relations T n = r~ ?1 T. Also the original time-dependent 
function, P(X, V, r), is replaced by an auxiliary function, P a (X, V, r , Ti, ..), depending on the r„, which are treated 
as independent variables. Once the equations corresponding to the various orders have been determined, one returns 
to the original time variable and to the original distribution. 

We begin by replacing the time derivative with respect to t by a sum of partial derivatives: 



(16) 



3 d 1 d 1 d 
The auxiliary function, P a {X, V, tq, t\, ...), is expanded as a series of powers of T -1 

oo ^ 

P a (X, V, T , n, t 2 , ...) = Y1 fI P a S) (*> V > to,ti,t 2 , ...). (17) 

s=0 

Each term Pa is then projected over the eigenfunctions H Ul 

oc 

p^(x, v, To ,n, ...) = 4 s) (x, T o, n,T 2) ...)H u (x, v). (is) 

i/=0 

At this stage one substitutes the time derivative l|16(l and expressions (|1TI) - (|18I) into Eq. l|10|l . and identifying terms 

1 .... . . (s) 

of the same order in T in the equations one obtains a hierarchy of relations between the amplitudes Vv • The 
advantage of the method over the naive perturbation theory is that secular divergences can be eliminated at each 
order of perturbation theory and thus uniform convergence is achieved. 

The most important function in all this expansion is . Its physical meaning becomes clear if one integrates 
over V the auxiliary density probability function P a (X, V, tq,ti,T2, —)■ Assuming, as in the standard expansion, 
that Vo"'' = for v > 1 one sees that ipo'' > is the spatial density probability function. We shall consider the partial 
derivative, d Ti ^\ and obtain the evolution equation for ■0q°' ) to order T~ l . For the sake of completeness we briefly 
show how the method works. To order T a one finds: 



£fp[5>1 0) ^] =0, (19) 



and concludes that only the amplitude is non-zero. 

We consider, next, terms of order T _1 and write (note that on the l.h.s we always have the application of the 
Fokker-Planck operator to the auxiliary function to the order of interest, in this case s = 1) : 

+ 24 1] H 2 + 34 1] H 3 + E,> 3 vi/)Ph u ) = H d To ^ 0) (20) 
+ ^T(x)H 1 [d x ^ ) - ^ + ^4 0) ] + i^H 3 4°\ 
By equating the coefficients multiplying Hq we find the relation: 

d 4 0) 



<9t 



= 0, (21) 



i.e. the amplitude ip^ is not a function of tq. We also obtain for the amplitude ip± the equations 



^^^ a ^ + 7mf''-M^- (22) 



The remaining amplitudes are 

1 T'(X) 



4 1} - 0, 4 1} = ^Tfjff o°\ 4 1] =0if » > 3- (23) 
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By iterating the method we obtain the correction of order T 1 to the evolution equation of j which reads 



and using Eq. (|22|l we can rewrite it as: 



d T1 4 0) = d x \d x (T^) - f4 0) 



(24) 



(25) 



Finally, we recover the original time variable, r, and obtain the most important result of this work, that is, the 
evolution equation for the spatial probability density: 



(o) _ 



a x (T^ 0) )-F^ 0) l+o(i/r 2 ) 



(26) 



With 0(l/r 2 ) we mean terms of or 1/T 2 and superior. The first comment that arises from Eq. is that it is 
of the Ito type if one uses this prescription for the overdamped Langevin equation. However, here we did not use 
any prescription for the Langevin equation, and Eq. I120H was derived from an expansion of the Kramers equation. 
Moreover, the generalization of Eq. i|26|) to arbitrary spatial dimensions is straightforward. 

Using a more usual notation for the spatial density probability function, and following the line of reasoning of 
Marconi and Tarazona™, one could write the following governing equation 



dp{x,t) 
dt 



= V 



[vD(x)p{x,t) 



1 , t yjv 5Fni\p] 
m-f ' Sp(x,t) 



V P . 



*(*)]}, 



(27) 



where D{x) — kBT(x)/m^ and F n i is the excess over the ideal gas contribution to the free energy^. The second term 
renormalizes the effective diffusion constant of the particles and may become important at higher concentrations. 

In the next section we proceed to check numerically our results and to enter in a possible discussion about the 
Ito-Stratonovich dilemma. 



IV. NUMERICAL RESULTS 



This section is devoted to numerically check our analytical results. The way to proceed is to compare the results 
of numerical simulations of a system of N particles following the dynamics given by Eqs. (JIJ or rather Eq. |JS}. We 
shall compare these with the analytical solution, if available, of the corresponding density equation: Eq. (|26[) for the 
first case, or the Smoluchowski in the Ito and Stratonovich prescriptions, Eqs. I|6I7|) . for the second case. 

In order to obtain a simple analytical solution of the density equation we assume a vanishing external force and 
molecular field, F — 0, particles have mass one, the system is one-dimensional and particles positions restrict to the 
interval [0, 1]. We also assume that the temperature field has a linear spatial dependence T(X) = Tq + T\X, with 
To and T± positive constants. Finally, we impose zero flux boundary conditions of the system at the extremes of the 
spatial interval, i.e., there are infinite walls at these extremes where the particles rebound and cannot scape from the 
system. 

The equation for the probability density that we have derived from the Kramers equation with the multiscale 
approach in an expansion in T, Eq. (|26[) . is now (for aesthetic reasons we denote in the following ip^ {X, r) = p(X, r)) 



d T p(X,r) = ±0* (T(A» = -d x J(X), (28) 

where J is the flux. The boundary conditions read J(0) = J(L) = 0. The stationary solution is given by: 

NTi 

pS{X) = (T + TxX) In [(T + TiZ)/T ] ' (29) 

where we have used the normalization condition dXp s {X) = N . 

On the other side, in the overdamped limit one has to distinguish between the two prescriptions. In this paper we 
discuss the two more standard: Ito and Stratonovich. The Ito-Smoluchowski equation is also given by Eq. I|28|) and 
the stationary solution is obviously given by Eq. (|29|l . A very important difference here is that it has not been derived 
as an expansion of inverse powers of T, and thus it is exact and no higher order corrections should appear. In the 
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FIG. 1: (color online) Plot of the stationary density of particles for different values of F. The particle dynamics correspond 
to the Langevin equation Eq. and the density of particles has been computed by making a coarse-graining in the particle 
positions. With the solid black line we plot the analytical solution Eg. 1291 that fits well for large F values. The other parameters 
used in the simulation are dt — 0.0001 and the cell size for the coarse-graining dx — 1/1000. 



Stratonovich prescription the density equation, which is also exact and not an approximation to order F , is given 
by: 

d T p{X,T) = IdxiT^dxT^p). (30) 
With the same boundary conditions the stationary Stratonovich solution is given by: 

NTl 

Pstrat ~ (T + T 1 X)i/2[2(VT + TiL - V%)] 

Now we show the numerical results for the Langevin particle dynamics. We begin with the underdamped dynamics, 
and consider 2000 particles evolving according to Eq. Q). We take To = 0-1 and 1\ = 10 and all our measurements 
are performed when the system is stationary. Average values are obtained by sampling 200 different realizations. A 
stationary density of particles is defined at every spatial point by counting the number of particles in every cell of size 
dx = 1/1000. The time step used is dt = 0.0001. For different values of T the results are shown in Fig.^ One can 
see that for large values of T the comparison with the analytical result, Eq. 1)29(1 . is excellent. Obviously, the smaller 
the value of T the larger the corrections that we are neglecting and the agreement worsens. 

In Fig. [5] we show the results obtained from the simulations of the overdamped dynamics, Eq. (J5J, where we have 
also coarse-grained the distribution of particles to compute a density of particles. Interpreting it like Ito the results 
are in perfect agreement with the Ito-Smoluchowski stationary solution, Ea. (|29|l . This is shown in the left panel of the 
figure for different values of T. When it is understood in the Stratonovich calculus, for which we have to add a new 
term T 2 !^ in the r.h.s. of Eq. (jjjj to make the numerical simulation, again the results coincide with the stationary 
density solution in the Stratonovich calculus, Eq. (|31|l . as shown in the right panel of Fig. [5] for various values V. 
A very important difference encountered in these results with respect to those in Fig. ^ is their independence with 
respect to the value of T: for any of its values the agreement with the stationary density is always excellent. In fact, 
one can realize that once the limit has been performed, the value of T only renormalises the time scale of the problem. 
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FIG. 2: (color online) Coarse-grained density field obtained from the overdamped Langevin dynamics Eq. [5] Left panel 
corresponds to the Ito interpretation of this equation, and right to the Stratonovich one. We take different values of the 
parameter V as indicated in the legend box. With the solid line we plot the stationary density solutions as given by Eqs. 1291 
and 1301 for Ito and Stratonovich respectively. The rest of the parameters used in the simulation are dt = 0.0001 and the cell 
size to make the coarse-graining of the position of the particles is dx — 1/1000. 



V. CONCLUSIONS 



In the present paper, we have derived the form of the governing equation for the probability density of interacting 
Brownian particles in the case of an inhomogeneous medium, whose temperature varies in space. The interparticle in- 
teractions have been treated in a mean-field framework. We have shown, using the Klein-Kramers equation associated 
to the fully inertial stochastic dynamics, that, to leading order, one obtains a Smoluchowski equation for the particle 
distribution which has the same form as the equation obtained by applying the Ito prescription to the overdamped 
dynamics. 
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